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1.0  INTRODUCTION 


The  formation  of  condensed  gases  on  cryogenically  cooled  surfaces 
has  recently  received  attention  in  regard  to  how  the  resulting  coating 
affects  the  reflection  and  transmission  of  visible  and  infra-red  (IR) 
radiation.  Sensor  testing  requires  knowledge  of  condensed  gas  coating 
reflectances  and  transmittances  because  film  deposits  formed  on  the 
sensor  system  transfer  optics  (which  are  cryogenically  cooled)  can  af¬ 
fect  their  absolute  throughput.  These  condensed  gas  deposits  both  ab¬ 
sorb  and  scatter  visible  and  IR  radiation.  Scattering  and  absorption 
are  functions  of  both  wavelength  and  coating  thickness.  In  addition,  if 
collimated  light  is  incident  upon  optical  surfaces  coated  with  condensed 
gases,  the  transmitted  and  reflected  light  will  no  longer  be  collimated. 
Thus  it  becomes  important  to  know  the  angular  distribution  of  the  re¬ 
flected  and  transmitted  energy.  Presently,  work  is  under  way  to  ex¬ 
perimentally  determine  the  effects  of  gas  deposits  formed  on  cryogeni¬ 
cally  cooled  (20°K)  lenses,  filters,  windows,  mirrors,  and  opaque  sub¬ 
strates  such  as  black  paint  and  stainless  steel. 

Experimental  data  are  important  in  determining  the  effects  of  coat¬ 
ing  thickness  and  viewing  angle,  as  well  as  wavelength,  upon  reflect¬ 
ance  and  transmittance.  However,  it  is  useful  to  have  analytical  or 
theoretical  expressions  which  mathematically  model  the  reflectance 
and  transmittance  behavior.  First,  the  analytical  expressions  are 
valuable  in  understanding  the  reasons  behind  the  behavior  of  the  ex¬ 
perimental  data.  (Also,  the  theoretical  expressions  can  be  used  to  de¬ 
termine  the  magnitude  of  the  various  contributing  processes  such  as 
scattering,  absorption,  substrate  effects,  front  interface  reflectance, 
etc.  In  this  regard,  the  theory  really  complements  the  experimental 
data.  )  Second,  when  the  theory  and  the  data  both  show  the  same  quan¬ 
titative  behavior  (as  in  Ref.  1),  then  the  theory  and  experimental  data 
may  be  used  to  determine  the  thickness  and  important  optical  properties 
of  the  condensed  gas  coating  such  as  refractive  index,  absorption  index, 
scattering  coefficient,  and  absorption  coefficient.  When  scattering  and 
absorption  are  both  present,  the  optical  properties  and  thickness  may 
be  obtained  through  using  experimental  data  in  conjunction  with  the  so¬ 
lution  to  the  radiative  transport  equation  (as  was  done  in  Ref.  1). 

The  results  obtained  in  Ref.  1  were  determined  by  solving  the  ra¬ 
diative  transport  equation  by  the  Milne  predictor- corrector  method. 

This  method  works,  but  it  is  cumbersome  to  use  and  requires  a  signif¬ 
icant  amount  of  computer  time.  In  an  attempt  to  find  a  simpler  and 
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faster  technique,  three  different  methods  for  solving  the  radiative 
transport  equation  are  presented  and  compared.  First,  the  transport 
equation  is  solved  iteratively  by  the  Milne  predictor- corrector  finite- 
difference  method.  Second,  an  integral  equation  formulation  for  the 
source  function  is  solved  by  successive  approximations;  knowing  the 
source  function  leads  directly  to  obtaining  the  intensity  as  a  function  of 
viewing  angle  and  optical  thickness.  Third,  the  Chandrasekhar  eigen¬ 
value,  eigenvector  solution  is  shown  in  comparison  with  the  other  two 
techniques.  It  is  important  to  have  several  solution  techniques  be¬ 
cause  each  technique  possesses  a  distinct  advantage  over  the  others. 

The  Chandrasekhar  method  is  simple  and  fast,  but  it  is  applicable  only 
when  the  absorption  and  scattering  coefficients  are  not  functions  of 
coating  thickness.  The  other  two  techniques  are  still  applicable  even 
when  the  optical  properties  vary  with  coating  thickness. 

In  choosing  three  techniques  for  comparison,  it  is  useful  to  review 
some  previous  work  requiring  solution  of  the  radiative  transport  equa¬ 
tion.  Wolf  (Ref.  2)  numerically  solved  the  transformed  (by  discrete 
ordinates)  transport  equation  via  Simpson's  rule.  After  reducing  the 
transport  equation  to  a  system  of  linear  ordinary  differential  equations. 
Wolf  obtained  solutions  for  problems  involving  absorbing,  emitting, 
and  scattering  media  with  an  arbitrary  specified  temperature  profile. 

In  Ref.  1  the  Milne  predictor-corrector  numerical  scheme  was  used 
to  investigate  the  reflectance  behavior  of  absorbing  and  internally  scat¬ 
tering  cryodeposits. 

The  formulation  of  an  expression  for  the  source  function  was  per¬ 
formed  by  Merriam  (Ref.  3).  Merriam's  technique  leads  to  an  integral 
equation,  the  solution  of  which  describes  the  source  function.  Once  the 
source  function  is  known,  the  radiant  intensity  behavior  can  be  directly 
computed. 

Besides  the  integral  equation  formulation  and  the  numerical  solu¬ 
tion,  other  solutions  to  the  transport  equation  have  been  obtained  through 
the  computation  of  eigenvalues  and  eigenvectors.  When  a  Gaussian  quad¬ 
rature  (discrete  ordinates)  is  used  to  approximate  the  integral  term  in 
the  transport  equation,  a  system  of  simultaneous  differential  equations 
is  obtained.  The  eigenvalues  and  eigenvectors  associated  with  the  coef¬ 
ficient  matrix  of  the  system  of  differential  equations  can  be  used  to  ob¬ 
tain  the  homogeneous  solutions;  emission  gives  rise  to  particular  solu¬ 
tions.  Once  the  general  solution  is  known,  the  boundary  conditions  are 
enforced  to  evaluate  the  integration  constants.  Hottel,  et  al.  (Ref.  4) 
employed  the  method  of  discrete  coordinates  to  compute  the  biangular 
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reflectance  from  an  absorbing  and  anisotropically  scattering  medium. 

The  system  of  differential  equations  was  solved  by  computing  the  eigen¬ 
values  and  eigenvectors.  Merriam  (Ref.  3)  has  also  made  use  of  the 
method  of  discrete  ordinates;  eigenvalues  and  eigenvectors  were  com¬ 
puted  by  the  "method  of  Danilevsky"  as  described  by  Faddeeva  (Ref.  5). 
Hsia  and  Love  (Ref.  6)  also  have  published  a  computational  method  of 
monochromatic  heat  transfer  in  the  plane,  one- dimensional  case  of  a 
parallel  atmosphere  separated  by  a  cloud  of  particles.  Applying  the 
method  of  discrete  ordinates,  they  solved  the  resulting  set  of  differen¬ 
tial  equations  by  obtaining  eigenvalues  and  eigenvectors  utilizing  the 
method  of  idempotents.  Another  particularly  attractive  method  for 
the  computation  of  eigenvalues  and  eigenvectors  for  solution  of  the 
transport  equation  in  association  with  both  isotropic  and  anisotropic 
scattering  is  that  of  Chandrasekhar  (Ref.  7).  Chandrasekhar's  method 
yields  an  equation  which  readily  permits  the  extraction  of  the  eigen¬ 
values  instead  of  requiring,  as  in  Ref.  5,  the  determination  of  the  char¬ 
acteristic  polynominal  before  attempting  computation  of  the  eigenvalues. 
In  addition  to  the  works  mentioned  above  there  are  numerous  other  solu¬ 
tion  methods  such  as  considering  radiant  transport  through  optically 
thick  media  to  be  a  diffusion  process  (Ref.  8).  Also,  the  transport  equa¬ 
tion  has  been  solved  by  Callis  via  the  method  of  characteristics  (Ref.  9). 

As  mentioned  earlier,  the  three  techniques  chosen  for  study  are 
the  Chandrasekhar  method,  the  source  function  integral  equation  formu¬ 
lation,  and  a  numerical  solution  by  the  Milne  predictor- corrector 
scheme.  These  three  were  chosen  because  they  are  easy  to  understand 
and  simple  to  apply,  and  because  they  serve  as  useful  tools  to  the  engi¬ 
neer.  Also,  it  is  useful  to  have  several  different  solution  techniques 
available,  as  one  technique  may  possess  a  distinct  advantage  over  the 
others.  In  order  to  form  a  basis  for  comparison  of  the  three  solution 
techniques,  a  standard  problem  has  been  selected  to  which  the  different 
methods  will  be  applied.  The  problem  chosen  is  that  of  an  absorbing 
and  isotropically  scattering  dielectric  medium  which  will  be  more  fully 
described  in  the  next  section  of  this  report.  Isotropic  scattering  was 
assumed  in  order  to  avoid  perplexion  with  a  highly  sophisticated  scat¬ 
tering  function.  However,  all  three  methods  presented  can  be  extended 
for  application  to  anisotropically  scattering  media.  Emission  was  con¬ 
sidered  negligible  in  order  not  to  require  particular  solutions  nor  to 
have  to  choose  a  specific  temperature  profile  or  temperature  level. 
Again,  the  three  techniques  discussed  here  can  easily  be  extended  to 
include  emission  once  the  temperature  profile  and  level  of  temperature 
have  been  prescribed. 
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The  three  solution  proceduresare  compared  with  one  another,  and 
the  Chandrasekhar  method  is  used  to  make  a  comparison  between  re¬ 
sults  obtained  by  single  and  double  Gaussian  quadratures.  Hottel,  et  al. 
(Ref.  4)  have  used  both  double  and  single  Gaussian  quadrature  approxi¬ 
mations  in  response  to  the  proposal  in  Ref.  10  that,  because  of  the  dis¬ 
continuity  of  the  intensity  for  p  =  0  at  the  bounding  interfaces,  the  quad¬ 
rature  formulas  should  be  applied  separately  in  the  p  ranges  of  -1  to 
0  and  0  to  1.  They  employed  as  many  as  20  double  Gaussian  quadrature 
directions  (points)  and  24  single  Gaussian  quadrature  directions  and  ob¬ 
served  no  appreciable  difference  between  the  two  types  of  quadratures. 
Hsia  and  Love  (Ref.  6)  have  used  the  double  Gaussian  quadrature  (8 
points)  in  their  work,  whereas  Wolf  (Ref.  2,  8  points),  Roux  (Ref.  1, 

10  points),  Merriam  (Ref.  3,  24  points),  and  Chandrasekhar  (Ref.  7) 
have  all  employed  the  single  Gaussian  quadrature.  One  purpose  of  this 
paper  is  to  determine  whether  there  is  any  computational  advantage 
(with  respect  to  accuracy,  convenience,  or  computer  time)  of  using  the 
single  as  compared  to  the  double  Gaussian  quadrature.  In  addition,  for 
each  quadrature  approximation  the  influence  of  the  number  of  quadra¬ 
ture  points  is  investigated;  as  many  as  96  points  were  used  for  the 
single  Gaussian  quadrature,  and  as  many  as  48  points  were  utilized 
for  the  double  Gaussian  quadrature. 


2.0  STATEMENT  OF  THE  PROBLEM 


In  order  to  compare  the  three  techniques  for  solving  the  transport 
equation  it  was  necessary  to  select  a  sample  problem  to  which  they 
would  be  applied.  The  geometry  and  coordinate  system  for  the  selected 
problem  are  shown  in  Fig.  1.  The  nomenclature  employed  in  Fig.  1  is 
essentially  the  same  as  that  used  in  Refs.  1  and  3;  regions  1,  2,  and 
3  respectively  represent  vacuum,  radiatively  participating  medium, 
and  opaque  substrate.  The  radiant  intensity  (IQ)  incident  on  the  radi¬ 
atively  participating  (dielectric)  medium  from  vacuum  is  taken  to  be  dif¬ 
fusely  distributed.  The  dielectric  coating  is  partially  transparent,  and 
the  opaque  substrate  is  considered  to  be  either  a  conductor  or  a  di¬ 
electric  having  negligible  emission.  The  vacuum- coating  interface  is 
assumed  to  reflect  radiant  intensity  in  accordance  with  Fresnel's  equa¬ 
tions  and  to  transmit  radiant  intensity  according  to  Snell's  law.  The 
substrate  is  taken  to  be  a  Fresnel  reflector.  Since  the  interface  re¬ 
flectance  for  radiant  intensity  traveling  from  vacuum  to  coating  and  the 
interface  reflectance  for  intensity  traversing  from  coating  to  vacuum  are 
not  the  same  for  equal  angles  of  incidence,  the  two  interface  reflectances 
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are  designated  as  pi2  and  P21*  respectively.  Also,  the  interface  re¬ 
flectance  for  intensity  incident  on  the  substrate  is  designated  as  P23. 
Further,  the  interior  of  the  radiatively  participating  coating  is  consid¬ 
ered  to  be  absorbing  and  isotropically  scattering  with  negligible  emis¬ 
sion.  Furthermore,  the  absorption  and  scattering  coefficients  (and 
hence  albedo)  are  considered  not  to  be  functions  of  the  local  position 
y  in  the  participating  film.  In  addition,  the  radiative  intensity  is  as¬ 
sumed  to  be  axisymmetric;  this  means  that  the  intensity  field  is  de¬ 
pendent  on  the  polar  angle  6  but  not  on  the  aximuthal  angle  <f>  .  The  sam 
pie  problem  thus  defined  via  Fig.  1  was  used  to  compare  the  three  so¬ 
lution  procedures  through  the  prediction  by  each  technique  of  the  mono¬ 
chromatic  hemispherical- directional  reflectance. 
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3.0  ANALYSIS 


Now  that  the  sample  problem  has  been  physically  defined,  the  math¬ 
ematical  presentation  of  the  governing  equations  will  be  given.  After 
the  transport  equation  and  boundary  conditions  are  shown,  each  of  the 
three  solution  techniques  will  be  applied  to  the  problem  as  a  subsection 
of  the  analysis.  The  radiative  transport  equation  subject  to  the  above 
stated  assumptions  is 

L.dHl  .1  (1) 

dr  n  -1 

where  i(T,ju)  is  the  local  radiant  intensity,  I(t,/u),  nondimensionalized 
by  the  diffusely  incident  intensity  I0  (Fig.  1);  /i  is  the  cosine  of  the  in¬ 
ternal  polar  angle  0  defining  the  direction  of  I(t,p);  W  =  cr/(a+k)  is  the 
albedo  parameter;  and  r  is  the  local  optical  depth,  which  is  related  to 
the  position  coordinate  y  by 


r  =  f\o  +  k)dy  (2) 

o 

with  a  and  k  being  the  scattering  and  absorption  coefficients,  respec¬ 
tively.  As  previously  stated,  the  coating  boundaries  were  considered 
to  have  only  regular  Fresnel  reflection  and  refraction;  expressed  math¬ 
ematically,  these  boundary  conditions  are 

=  p2](/z)  i(r0,^)  -  jl  -  ^ 

for  the  vacuum- coating  interface  and 

i(0,p)  =  P23W  i(0,-jx)  (4) 

for  the  coating- substrate  interface  where  n2  is  the  refractive  index  of 
the  coating  and  pt1  is  the  cos  0j,  with  0^  being  the  polar  angle  external 
to  the  radiatively  participating  medium  (Fig.  1).  The  directions  ju  and 
/u1  are  related  by  Snell's  law, 

p1  =  [l  — (l  — /i2)n|f  (5) 

and  t0  =  r(y=L)  where  L  is  the  coating  geometrical  thickness, 

3.1  METHOD  OF  CHANDRASEKHAR 

The  basic  objective  of  the  method  of  Chandrasekhar  is  computation 
of  the  eigenvalues  and  eigenvectors  of  the  coefficient  matrix  associated 
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with  the  system  of  simultaneous  ordinary  linear  differential  equations 
resulting  from  the  use  of  discrete  ordinates.  Once  the  eigenvalues  and 
eigenvectors  are  determined,  the  homogeneous  solution  is  known;  the 
integration  constants  are  then  directly  obtainable  from  the  boundary 
conditions.  Before  the  eigenvalues  and  eigenvectors  can  be  computed, 
the  transport  equation,  Eq.  (1),  together  with  the  boundary  conditions, 
Eqs.  (3)  and  (4),  must  be  first  transformed  into  a  system  of  differential 
equations  by  the  method  of  discrete  ordinates.  This  consists  of  replac¬ 
ing  the  integral  term  in  Eq.  (1)  by  a  Gaussian  quadrature  of  the  form 


f{x)dx 


-  £  ai  f(xi> 


(6) 


where  Xj  is  the  quadrature  points  (discrete  ordinates),  aj  is  the  quadra¬ 
ture  weights,  and  p  {which  is  an  even  integer)  is  the  order  of  the  quad¬ 
rature  approximation.  Replacing  the  integral  term  in  Eq.  (1)  by  Eq. 

(6)  yields  a  system  of  p  simultaneous  differential  equations 


di(r,fi|0 

dr 


+  — —  £  a-i(r,/ij),  £-1 . p 

2/ig  j=l  >  ^ 


with  boundary  conditions 


i(r0.-fzg)  =  p2 iWWVMe)  +  [l -P]2(M£)n2- E=1 . p/2 


and 


(7) 


(8) 


i(0,j/g)  =  *  G=l,...,p/2 


(9) 


The  coefficient  matrix  of  Eq.  (7)  can  be  written  as 

B£j  =  (Waj/2  —  Sgp/ up  ,  E,j  =  l,...,p  (10) 

where  6jg j  is  the  Kronecker  delta.  This  matrix  has  exactly  the  same 
form  whether  the  single  Gaussian  quadrature  is  employed  over  the  in¬ 
terval  -1  <_1  as  shown  in  Eq.  (6)  or  whether  the  double  Gaussian 

quadrature  is  used  separately  in  the  regions  0  to  1  and  -1  to  0.  For 
the  latter  case,  the  double  Gaussian  quadrature  of  order  p/2  is  used 
twice,  thus  giving  a~p  x  p  coefficient  matrix  of  the  form  in  Eq.  (10). 
The  quadrature  points  and  weights  are  different  for  the  two  cases,  but 
in  both  instances  they  possess  the  common  properties 

J*k  =  “Fp+i—k  ’  ak  =  ap+l_k  ’  k  =  l,...,p/2  (11) 
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It  has  been  shown  {Ref.  11)  that  the  p  eigenvalues  associated  with 
the  p  x  p  coefficient  matrix,  Eq.  (10),  are  the  p  values  of  X  which  sat¬ 
isfy  the  expression 


(12) 


Note  that  Eq.  (12)  is  a  function  of  X^.  It  is  known  that  the  roots  X^  are 
positive  or  zero  as  shown  in  Fig.  2.  Thus,  the  eigenvalues  occur  as 
plus  and  minus  pairs.  When  W  =  0,  the  first  factor  of  Eq.  (12)  must 
be  zero,  giving  the  p  eigenvalues  as  Xj  =  ±1  //ij,  j  =  1,  .  .  . ,  p/2.  If 


1  1  1 

A  A  A 

Note.-  Plot  is  not  to  scale. 

Figure  2.  Sketch  illustrating  graphical  solution  for  the 
eigenvalues,  p  =  6,  Ref.  1 1 . 
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W  f  0,  the  second  factor  must  be  zero.  Setting  the  second  factor  equal 
to  zero  and  simplifying,  one  finds  the  eigenvalues  must  satisfy  the 
equation 


gU2) 


p/2  BjM?X2 


which  is  equivalent  to  that  obtained  in  Ref.  7, 


(13) 


p/2  a. 


I 

w 


(14) 


if  one  notes  that  the  term  - 1  on  the  right  side  of  Eq.  (13)  is  equivalent 

p/2 

to  -  £  a-i.  When  W  =  1.0,  the  right  side  of  Eq.  (13)  is  zero,  and  it  is  seen 

j=l  J 

that  X2  =  0  is  a  root.  Thus,  zero  is  an  eigenvalue  with  multiplicity  of 
two  and  requires  a  special  form  of  the  homogeneous  solution  in  order 
to  insure  linearly  independent  solutions.  Plotting  simultaneously  g(X2) 
and  1/W-l  versus  X2  yields  the  solutions  X2  where  g(X2)  and  (1/W-l)  in¬ 
tersect.  Figure  2  is  an  illustrative  plot  for  p  =  6,  but  the  analysis  is 
valid  for  any  even  integer  p.  It  is  seen  from  Fig.  2  that  the  values  X2 
satisfying  Eq.  (13)  must  be  positive  or  zero.  Thus,  all  the  eigenvalues 
must  be  real  and  must  occur  in  plus  and  minus  pairs  except  for  W  =  1, 
where  two  of  the  eigenvalues  are  equal  to  zero.  Also,  from  Fig.  2  one 
can  readily  establish  the  bounds  for  the  roots  X2.  For  this  example, 

0  <  X|  £  1/p^  <  x|  <L  1/^1  <  ^3  1/m §•  Thus  all  the  roots  X2  have  indi¬ 

vidual  bounds.  Using  these  bounds,  the  numerical  solution  for  the  roots 
of  Eq.  (13)  may  be  easily  obtained.  The  homogeneous  solution  of  Eq. 

(7)  for  W  ^  1.0  is  given  by 


i(r,p£) 


P/2 
=  1 
i=i 


i  -x:1 


1 "  fyt 


[cj(l  -  Xyxf;)eAjr 


cp+l_j(1  -  Xj/i£)e_X3 


'] 


1 >  ■  *.»F 


(15) 


where  the  c’s  are  the  p  number  of  integration  constants  determined  from 
the  boundary  conditions.  Substitution  of  Eq.  (15)  into  Eqs.  (8)  and  (9) 
yields  a  system  of  p  nonhomogeneous  linear  algebraic  equations  to  be 
solved  for  the  p  values  of  c;  use  of  the  Gauss-Jordon  method  or  Cholesky 
method  (see  Ref.  12)  readily  allows  determination  of  the  integration 
constants.  For  the  special  case  of  W  =  1.0,  the  double  eigenvalue  X^  = 
Xp=0  causes  Eq.  (15)  to  have  the  special  form  (Ref.  13) 

(16) 


i(r,pg)  = 


+  c 


A  -  w* 


P/2 

j=2  j  _ 


-Y 


Cp+H(1 


— -A:r" 

Yf)e  J 


13 


A  EDC-TR -73-200 


Chandrasekhar  (Ref.  7)  also  has  expressions  equivalent  to  Eqs. 
and  (16)  if  one  defines  the  Chandrasekhar  integration  constants, 

c£+i-j'  as 


(15) 


ci 


and 


cj  = 


and 


CP+l-i  “  (1  'WVl-j 


(17) 


It  should  be  noted  that  Eqs.  (15)  and  (16)  yield  the  solution  of  Eq.  (7) 
in  the  p  quadrature  directions.  If  information  is  required  in  directions 
other  than  the  quadrature  directions,  then  interpolation  must  be  used 
among  the  solutions  in  the  quadrature  directions.  A  final  comment 
stressing  the  simplicity  of  this  technique  is  that  once  the  eigenvalues 
are  obtained  from  Eq.  (13)  (which  is  easily  done  since  their  bounds  are 
known),  the  homogeneous  solution  is  immediately  given  by  Eq.  (15)  or, 
if  W  =  1.0,  by  Eq.  (16).  Only  the  integration  constants  remain  to  be 
found,  and  these  can  readily  be  obtained  by  standard  computer  library 
routines  for  solving  systems  of  nonhomogeneous  linear  algebraic 
equations. 

The  eigenvectors  have  already  been  included  in  Eqs.  (15)  and  (16), 
and  they  are  functions  only  of  the  quadrature  points  and  the  eigenvalues. 
Information  on  the  derivation  of  the  eigenvectors  is  available  in  Refs. 

7,  11,  and  13. 


3.2  SOURCE  FUNCTION  FORMULATION 

Solution  of  the  transport  equation,  Eq.  (1),  along  with  the  boundary 
conditions,  Eqs.  (3)  and  (4),  through  the  source  function  formulation 
was  accomplished  by  solving  an  integral  equation  describing  the  source 
function.  The  intensity  field  was  then  computed  through  its  direct  re¬ 
lation  to  the  source  function.  In  order  to  arrive  at  the  integral  equation, 
let  Eq.  (1)  be  rewritten  as 


dr  n  ~  n 


where  4>(r)  is  the  normalized  source  function 


<t(r) 


W 


2 


0 


f  i(r,fi')d/i 


* 


(18) 


(19) 
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Solving  Eq.  (18)  for  the  normalized  intensity  yields  solutions  for  the 
intensity  in  the  upward  (positive  iu)  and  downward  (negative  n)  directions 
as 


♦  ft f 

o  r 


(20) 


and 


i =  i(r0,-fx)e  ^  h  /  —  (21) 

r  V 

where  now  in  Eqs.  (20)  and  (21)  0  <  ju  <  1.0.  From  Eqs.  120)  and  (21) 
it  is  observed  that  the  intensity  in  the  upward  and  downward  directions 
is  directly  related  to  the  normalized  source  function.  Hence,  once  the 
source  function  is  known,  the  intensity  field  can  be  directly  evaluated. 
The  procedure  to  be  followed  here  in  determining  the  source  function 
is  first  to  find  i(0,fi)  and  i(r0,  ~iu)  in  Eqs.  (20)  and  (21)  through  the 
boundary  conditions.  Once  i(0, jj)  and  i(TQ,  -ju)  are  known,  then  Eqs. 
(20)  and  (21)  are  substituted  into  Eq.  (19),  resulting  in  one  integral 
equation  to  be  solved  for  the  one  unknown,  $»(t). 

For  the  problem  defined  in  this  paper,  the  integral  equation  is  ob¬ 
tained  by  proceeding  as  outlined  above.  From  Eqs.  (20)  and  (21)  it  is 
found  that 


and 


i(r0,/i)  =  i(0,p)e 


TJV-  /”  ,  -fro-Wf*  cl  i 

+  f  <P(r)e  — 

6  M 


*(0 , — jx) 


(22) 


(23) 


Substituting  the  expressions  i (t0,/h)  and  i(0,  -ju)  from  Eqs.  (22)  and  (23) 
into  Eqs.  (3)  and  (4)  readily  yields  the  values  i(0,ju)  and  i(r0,  -fu)  needed 
in  Eqs.  (20)  and  (21).  With  i(0,/u)  and  i(T0, -ju)  known,  Eqs.  (20)  and 
(21)  are  substituted  into  Eq.  (19),  which,  after  the  indicated  algebra  is 
performed,  becomes  the  integral  equation  describing  4>(t), 
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m  =  r 


D<Vfd 


d H  +  f 


D(T0,p) 


+  /  °4>(t)  dt  +  /  °<D(t)  _ ! _ dt 


D(r0*t)  n 


-y2r  +r-i)/(i 

♦  I  °m  jg! - i  dt  +  ;  "®(t)  r 

o  o  D  (r0  ,/x)  poo 


— (2r  — r— t)/u 

r°  -  lP21^)e  d, 


-  ^  dt 

D(r0,^)  n 


r  -(2r -r+t)//i 

+  /  °  O(t-)  ^  dt 

°  °  D(ro,/d  P 


(24) 


where 

D(vp)  =  l-P21^P23^e"2r°/'1  <25> 

The  solution  of  Eq.  (24)  was  obtained  through  the  method  of  successive 
approximations.  At  first  glance  Eq.  (24)  appears  to  be  quite  compli¬ 
cated;  however,  closer  observation  reveals  that  this  is  not  so.  First 
notice  that  the  integrands  of  the  first  two  integrals  of  Eq.  (24)  are  all 
known  functions,  hence  these  integrals  only  need  to  be  evaluated  once. 
Since  their  values  can  readily  be  computed,  these  terms  do  not  have  to 
be  reevaluated  during  the  process  of  successive  approximations.  A 
similar  observation  is  noticeable  concerning  the  integration  with  re¬ 
spect  to  p.  in  .the  last  five  terms  on  the  right  side  of  Eq.  (24).  In  the 
last  five  terms  of  Eq.  (24)  all  the  functions  which  are  integrated  with 
respect  to  p  are  known  functions;  thus  the  integration  with  respect  to 
p  need  only  be  performed  once  during  the  process  of  successive  approxi¬ 
mations.  These  observations  can  result  in  a  significant  saving  of  com¬ 
puter  time,  but  a  large  amount  of  storage  is  required  since  the  values 
of  the  integrals  with  respect  to  p  in  the  last  five  terms  of  Eq.  (24)  must 
be  stored  at  each  t  and  t  combination  (t  occurs  in  the  exponent  of  the 
exponential  inside  the  integration  with  respect  to  p). 

The  solution  of  Eq.  (24)  was  determined  at  discrete  r  values  by  the 
method  of  successive  approximations.  The  integration  with  respect  to 
p  was  performed  by  using  10  Gaussian  quadrature  points.  The 
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integration  with  respect  to  t  was  accomplished  via  Simpson's  rule  with 
the  optical  thickness,  t0,  divided  into  50  equal  intervals.  The  initial 
guess  at  the  solution  for  $(r)  was  <J>(r )  =  0.0.  The  initial  guess  at  ♦,(t) 
is  substituted  into  the  right-hand  side  of  Eq.  (24),  and  an  improved 
estimate  of  <D(t)  is  computed.  This  process  is  repeated  until  converg¬ 
ence  is  considered  to  have  been  obtained.  In  this  work  the  successive 
approximations  were  repeated  until  the  maximum  difference  between 
two  successive  calculations  of  <E(r)  at  any  t  station  was  less  than  0.0001. 
Once  the  process  of  successive  approximations  has  converged,  the 
value  i(r0,/i)  can  be  directly  computed  from 

-  <26» 

The  integrals  in  Eq.  (26)  were  again  evaluated  via  Simpson's  rule.  The 
expression  i(T0, n),  as  will  be  shown,  is  needed  to  compute  the  hemi¬ 
spherical-directional  reflectance.  Note  that  Eq.  (26)  yields  informa¬ 
tion  in  the  quadrature  directions  and  in  all  other  directions;  no  inter¬ 
polation  is  required  to  obtain  information  in  the  nonquadrature  direc¬ 
tions,  as  is  necessary  for  the  Chandrasekhar  method  and  the  predictor- 
corrector  numerical  solution. 


3.3  NUMERICAL  SOLUTION 

After  the  transport  equation  and  boundary  conditions,  Eqs.  (1),  (3), 
and  (4),  were  transformed  by  the  use  of  discrete  ordinates  into  the  sys¬ 
tem  of  equations  given  in  Eqs.  (7)  through  (9),  solution  was  obtained  by 
the  Milne  predictor- corrector  method.  The  system  of  equations  to  be 
solved  is  the  same  system  as  that  solved  by  the  Chandrasekhar  method. 
The  primary  advantage  of  the  numerical  solution  (and  source  function 
formulation)  is  that  the  system  of  differential  equations  can  be  allowed 
to  have  variable  coefficients;  the  Chandrasekhar  method,  although 
powerful,  is  applicable  only  when  the  system  of  differential  equations 
has  constant  coefficients.  In  the  application  of  the  predictor- corrector 
method  to  the  system  of  equations,  the  number  of  steps  utilized  in  the 
numerical  calculations  was  varied  from  50  to  200.  The  mathematical 
expressions  for  the  predictor- corrector  method  adapted  to  i(r, n%), 

S.  -  1,  .  .  . ,  p/2  are 
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Predictor 


'(rn+|'W'  =  '(rn— 3  ’/*?)  +  Y 

"2  Y(rn’^  "  +  2  <rn-2'/^] 

(27) 

Corrector 

“1 

i(rn+r^  =  i(rn-I’W)  +  +  4  £  <  VW>)  +  Tr  ('n+r/^J 

(28) 

and  for  i(r .Mjg),  &  -  p/2  +  1,  .  . 

,  .  ,  p  are 

Predictor 

4h 

d  i  (I  i  d  i  | 

+  2  dr  rnH-2,/z?J 

(29) 

Corrector 

£  (fm+lrf  +  4Y(r"’’^)  +  Tt 

(30) 

where  n  =  4,  5,  .  . .  ,  N  and  m  - 

N-4,  N-3,  .  .  .  1  w'ith  K  being  the  number 

of  stations  employed  in  marching  from  0  to  rQ,  and  h  the  step  size, 
h  =  r0/(N-l). 

In  order  to  use  Eq.  (27),  one  must  know  the  derivatives  at  Tn,  Tn_1, 
and  Tn_2-  These  values  can  be  determined  directly  from  the  differential 
equations  only  after  the  corresponding  values  at  Tj,  T2,  T3,  and  t4  are 
known.  Hence,  a  starting  equation  such  as  the  modified  Euler  or  Runge- 
Kutta  technique  must  be  employed.  Using  the  modified  Euler  technique 
and  the  predictor- corrector  method,  the  system  of  differential  equa¬ 
tions,  Eq.  (7),  was  solved  iteratively.  Since  this  problem  is- boundary 
valued,  it  was  necessary  to  make  an  initial  guess  of  the  intensity  field 
at  the  top  (t=t0)  and  bottom  (t  =  0)  interfaces  before  a  starting  equation 
could  be  employed.  The  initial  guess  was  made  by  use  of  a  generalized 
version  of  the  dual  beam  approximation  discussed  in  Ref.  8.  This  gen¬ 
eralized  dual  beam  approximation  (Ref.  14)  not  only  allows  guessing 
the  intensity  field  at  t  =  0  and  t=t0  but  also  permits  guesses  at  any  of  the 
N  stations  associated  with  r  in  each  of  the  p  quadrature  directions  cor¬ 
responding  to  p.  The  quadrature  directions  of  the  intensity  traveling 
upward  toward  the  vacuum-deposit  interface  are  n j2  =  1,  ....  p/2,  and 
those  associated  with  intensity  traveling  toward  the  deposit-substrate 
interface  are  Pjg,  i=p/2+l,  . . . ,  p. 

The  actual  iteration  routine  used  was  a  forward  and  backward  in¬ 
tegration  scheme  which  proceeds  as  follows:  (1)  The  input  parameters 
such  as  refractive  indices,  optical  thickness,  and  albedo  are  chosen; 

(2)  the  interval  between  0  and  tq  is  divided  into  equal  intervals  at  N 
stations,  and  the  order,  p,  of  the  single  Gaussian  quadrature  is  selected; 

(3)  the  generalized  dual  beam  approximation  was  used  to  obtain  an  initial 
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guess  of  the  dimensionless  radiation  field  in  the  p  quadrature  directions 
at  every  t  station.  For  example,  t=t0,  an  initial  guess  was  obtained 
for  i(T0,ju^),  i =p/ 2+1,  .  .  . ,  p/2;  (4)  Eq.  (7),  4=p/ 2+1, .  .  . ,  p,  was  inte¬ 
grated  backward  (downward)  over  t  from  t=tq  to  t=0  by  holding  the 
guessed  values  of  i(r,  p_g),  i=T,  .  .  . ,  p / 2  constant  and  calculating  new 
values  of  i(T,Pjg),  i  =  p/  2+1,  .  .  . ,  p  by  use  of  Eqs.  (29)  and  (30);  (5)  after 
arriving  at  t  =  0,  the  boundary  condition;  Eq.  (9),  was  used  to  compute 
new  values  of  i(0,  p^j),  i=l,  ....  p/2;  (6)  Eq.  '(7),  i  =  1,  ...,p/2,  was 
integrated  forward  (upward)  from  t  =  0  to  t=t0  by  holding  the  newly  cal¬ 
culated  i(-r,p,g),  jtf=p/2+l,  .  . . ,  p  values  constant  and  computing  new  val¬ 
ues  of  i(T, n$),  JJ=1,  .  .  .,  p/2  through  the  use  of  Eqs.  (27)  and  (28)  (After 
arriving  at  t=tq  new  values  of  i(T0 ,p,g),  j£  =  1,  ....  p/2  were  available. ); 
(7)  the  values  of  IU^q.m^)!  ~  i(T0, M )0 1 jg * •  J0=1,  ■  •  • »  p/2  were  compared 
with  a  given  tolerance.  (In  this  case  the  maximum  value  had  to  be  less 
than  0.  001.  If  the  difference  in  any  quadrature  direction  was  greater 
than  the  given  tolerance,  the  boundary  condition,  Eq.  (8),  was  used  to 
find  new  values  for  i(+0,Mj g)>  J0=p/2+l,  ....  p,  and  the  entire  backward 
and  forward  marching  process  was  begun  again  at  step  4.  );  (8)  the 
procedure  in  steps  4  through  7  was  continuously  repeated  until  the  maxi¬ 
mum  difference  between  the  values  of  i(T0,p,g),  i  =  l,  .  .  . ,  p/2,  for  two 
successive  iterations  was  within  the  given  tolerance. 


3.4  DEFINITION  OF  phd{^) 


Before  proceeding  into  a  discussion  of  the  results,  it  is  necessary 
to  define  the  hemispherical- directional  reflectance.  The  hemispherical- 
directional  reflectance  is  the  ratio  of  the  intensity  reflected  from  an  in¬ 
finitesimal  area,  dA,  collected  in  a  specific  angular  direction  to  the  dif¬ 
fusely  incident  intensity  which  is  hemispherically  distributed.  Mathe¬ 
matically,  the  hemispherical- directional  reflectance  for  the  problem 
defined  in  Fig.  1  is 


Phdfo1)  =  + 


Kr0,fi) 


[1  -  p2](.p)] 


(31) 


2  1/9  i 

where  (1- l/n^)'  <  pi  <  1  and  p  is  related  top1  by  Eq.  (5).  Equation 

(31)  can  be  used  directly  in  conjunction  with  Eq.  (26)  to  determine  P^O^) 
through  the  source  function  formulation.  For  the  Chandrasekhar  tech¬ 
nique  and  the  predictor- corrector  method,  solution  to  the  transport 
equation  was  obtained  in  terms  of  the  dimensionless  intensity  in  the  p^ 
quadrature  directions.  The  hemispherical- directional  reflectance  in 
terms  of  the  quadrature  directions  becomes 


+  [l-p21W>J - o 

n 


(32) 
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where  (l-l/n^)1^2  <  —  1-»  -£  =  1»  •  •  • ,  p/  2  with  /uj  related  to  p  ^  by 

Snell’s  law,  Eq.  (5).  The  value  of  K^,/^)  in  Eq.  (32)  was  obtained 
either  from  the  predictor- corrector  solution  or  from  the  Chandrasekhar 
solution  via  Eq.  (15)  for  W  ^  1.0  or  via  Eq.  (16)  for  W  =  1.0. 


4.0  RESULTS 


The  results  obtained  from  the  three  solution  techniques  will  now  be 
compared.  Because  the  numerical  agreement  of  some  of  the  results 
was  so  good,  it  was  decided  that  the  results  should  be  illustrated  through 
the  use  of  tables  instead  of  by  showing  several  coincident  curves  on  a 
figure.  Tables  1  and  2  show  a  comparison  of  the  hemispherical- 
directional  reflectance  results  predicted  by  the  various  methods  for 
n2  =  1.2  and  tq  =  0.5  and  5.0  and  corresponding  to  the  various  /uj  quad¬ 
rature  directions.  The  results  shown  were  obtained  using  a  10-point 
single  Gaussian  quadrature.  Only  Pfod^^  in  three  of  the  five  upward 
directions  is  shown  since  the  other  two  directions  correspond  to  inten¬ 
sity  trapped  in  the  coating  because  of  being  incident  at  the  coating- 
vacuum  interface  at  angles  greater  than  the  critical  angle.  Also  shown 
in  Tables  1  and  2  are  the  results  obtained  by  the  method  of  Danilevsky 
(Ref.  5).  This  method  should  be  considered  as  a  "library”  program 
since  the  actual  programming  of  the  Danilevsky  technique  was  not  per¬ 
formed  by  the  authors  of  this  report. 


Table  1.  Comparison  of  Phdtju1 )  Results  for  the  Three  Solution 
Techniques:  n2  =  1.2,  W  =  0.4 


Substrate 

T0 

Chandrasekhar  Danilevsky 

Numerical 

Solution 

Source 

Function 

Black  Paint 

d  96220 

0  5 

d  04557 

0. 04557 

004657 

0  04774 

d 79850 

d  04998 

0.04998 

005008 

0.04895 

iyl.  48-id  00 

a  47403 

0.08522 

0.08522 

0  08795 

0.07732 

0.96220 

5.0 

0.06271 

0.06277 

0.06271 

007019 

0.79850 

0  06741 

0  06741 

006741 

007161 

d  47403 

010120 

0  10120 

0  10130 

0.09893 

Steel 

d  96220 

0  5 

02527 

0.2528 

02520 

0  2567 

a  79850 

0.2352 

02353 

0.2363 

02348 

n3  -  2. 48  -  i3. 43 

0.47403 

0.2202 

02202 

0.2272 

02107 

0.96220 

5.0 

006275 

006275 

006271 

007022 

d 79850 

0  06743 

006743 

006741 

007162 

d  47403 

0.10120 

010120 

0  10130 

009895 
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Table  1.  Concluded. 


Substrate 

P1 

T 

O 

Chandrasekhar  Danilevsky 

Numerical 

Solution 

Source 

Function 

Aluminum 

0. 96220 

as 

0.3663 

03662 

03676 

03721 

a  79850 

a  3374 

0.3374 

0  3373 

03380 

n3  •  L  44  -  i5. 32 

Q.  47403 

0.2965 

02965 

02974 

02865 

a  96220 

5.0 

0  06277 

006271 

006271 

007024 

0. 79850 

006744 

006744 

006741 

007164 

a  47403 

010130 

010130 

010130 

009896 

Copper 

a  96220 

0L5 

04400 

04397 

04418 

04475 

a.  79850 

04045 

0  4044 

04044 

04063 

n3*  082  -113.00 

a  47403 

03491 

0  3490 

03498 

0  3394 

0.96220 

5.0 

006278 

006278 

005271 

007025 

0.79850 

006744 

006744 

006741 

0.07165 

0. 47403 

010130 

0  10130 

010130 

0  09897 

Table  2.  Comparison  of  PhdfM1 )  Results  for  the  Three  Solution 
Techniques:  n2  =  1.2,  W  =  1.0 


„ 1 

Numerical 

Source 

Substrate 

“i 

To 

Chandrasekhar  Danilevsky 

Solution 

Function 

Black  Paint 

096220 

05 

0.1591 

01591 

01591 

01688 

079850 

01742 

01742 

01742 

0. 1698 

L48  -  iOOO 

0  47403 

0.2293 

02294 

02294 

01951 

096220 

5.0 

0  6903 

06903 

0  6861 

NC3 

079850 

07058 

07068 

07018 

047403 

07391 

0  7391 

0  7356 

Steel 

0  96220 

05 

05011 

05013 

05015 

0.5338 

079850 

05025 

05026 

05029 

05114 

n3  -  2.48  - 13.43 

047403 

05188 

05190 

0  5192 

04795 

096220 

5.0 

07454 

07455 

0  7407 

NC3 

0  79850 

0  7581 

07581 

0  7535 

047403 

0  7854 

07855 

07826 

Aluminum 

096220 

05 

0  7468 

0  7468 

07472 

0  7881 

0  79850 

0  7444 

07444 

07448 

0  7626 

rij  -  L  44  - 15. 32 

0  47403 

0  7453 

07453 

0  7457 

0  7243 

0  96220 

5.0 

08222 

08222 

0.8346 

NCa 

079850 

08310 

0  8310 

08428 

047403 

08501 

0.8501 

08606 

Copper 

096220 

0.5 

09615 

09609 

09619 

09951 

0  79850 

09606 

09603 

0  9611 

09727 

n3  ■  0  82  - 1 13, 00 

047403 

0  9596 

09593 

0.9601 

09460 

096220 

5.0 

09600 

0.9597 

09709 

NC3 

079850 

09620 

0.9617 

0  9723 

047403 

09663 

0.9660 

09754 

aNo  convergence  after  50  iterations  using  50,  100,  and  200  integration  steps  for  t0  end  10 
quadrature  points  for  the  Integration  with  respect  to  p  In  the  interval  0  <  p  <  1,  Eq.  (24). 
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The  Danilevsky  method  was  used  to  compute  the  coefficients  of  the 
characteristic  polynominal  from  the  coefficient  matrix,  Eq.  (10).  The 
eigenvalues  are  the  roots  of  the  characteristic  polynominal  and  were 
determined  by  the  Newton- Raphson  method;  then  the  eigenvectors  were 
constructed.  With  the  eigenvalues  and  eigenvectors  known,  the  inte¬ 
gration  constants  were  determined  via  the  Cholesky  method.  The 
Danilevsky  method,  as  used  here,  was  found  not  to  be  capable  of  hand¬ 
ling  a  coefficient  matrix  larger  than  16  x  16  since  the  trace  of  Eq.  (10) 
yielded  a  contradiction.  The  trace,  which  is  the  sum  of  the  roots 
(eigenvalues)  of  the  characteristic  polynominal,  was  no  longer  equal  to 
the  sum  of  the  diagonal  elements  of  the  coefficient  matrix,  Eq.  (10). 

The  results  obtained  by  the  Danilevsky  method  have  also  been  included 
in  Tables  1  and  2  as  an  additional  indication  of  the  accuracy  of  the  other 
results  shown. 

Table  1  corresponds  to  W  =  0.4,  and  Table  2  corresponds  to  results 
obtained  for  W  =  1.0.  The  agreement  of  the  results  predicted  by  the 
various  techniques,  as  illustrated  in  Table  1,  is  seen  to  be  good;  the 
agreement  obtained  by  the  Chandrasekhar  solution,  Danilevsky  method, 
and  numerical  method  is  excellent.  Table  2  (W  =  1.0)  also  shows  ex¬ 
cellent  agreement  among  the  Chandrasekhar  solution,  Danilevsky  method, 
and  numerical  solution;  however,  the  successive  approximations  solu¬ 
tion  to  the  integral  equation  did  not  converge  for  large  tq  values.  Table 
2  also  shows  poor  agreement  for  the  integral  equation  solution  at  small 
r0  values. 

The  Chandrasekhar  and  Danilevsky  methods,  as  expected,  were 
found  to  be  several  orders  of  magnitude  faster  than  the  predictor- 
corrector  method  or  integral  equation  solution.  The  Chandrasekhar 
and  Danilevsky  methods  have  the  same  common  disadvantage  in  that 
both  are  applicable  only  when  W  is  not  a  function  of  t.  The  Chandrase¬ 
khar  method  is  the  superior  of  the  two  because  it  permits  easier  and 
more  accurate  determination  of  the  eigenvalues,  it  is  faster  with  re¬ 
spect  to  computer  time,  and  (as  will  be  shown)  it  permits  much  higher 
order  quadrature  approximations  (96  x  96)  with  essentially  no  sacrifice 
of  accuracy.  The  Chandrasekhar  solution  (Eq.  15  or  16)  also  has  the 
eigenvectors  explicitly  known,  whereas  for  the  Danilevsky  method  the 
eigenvectors  must  be  repeatedly  constructed.  The  computation  of  the 
eigenvalues  by  the  Chandrasekhar  method  is  extremely  fast,  and  the 
PhdO-**)  is  directly  given  by  Eqs.  (15)  or  (16)  and  (32). 

For  W  =  0.4,  both  the  predictor-corrector  method  and  the  succes¬ 
sive  approximations  technique  (for  the  integral  equation)  were  found  to 
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converge  very  rapidly;  the  amount  of  computer  time  required  to  obtain 
a  solution  was  about  the  same  for  both  cases,  with  the  predictor- 
corrector  solution  being  slightly  faster.  The  predictor-corrector 
method  was  probably  faster  because  the  dual  beam  approximation  (Ref. 
14)  was  employed  as  an  initial  guess  criterion.  For  all  practical  pur¬ 
poses,  however,  the  two  methods  are  equivalent  with  respect  to  com¬ 
puter  time.  The  successive  approximations  procedure  required  about 
6  iterations  to  converge,  whereas  the  numerical  solution  required  less 
than  6  and  often  only  3  iterations.  The  predictor-corrector  solution 
was  found  to  have  better  agreement  with  the  Chandrasekhar  results. 

For  W  =  1.0,  the  Chandrasekhar  and  Danilevsky  methods  were  again 
very  fast.  However,  the  predictor-corrector  technique  was  found  now 
to  require  about  12  iterations  before  converging.  This  was  probably 
because  the  initial  guess  method  did  not  yield  as  good  an  initial  guess 
at  W  =  1.0  as  at  W  =  0.4.  In  spite  of  the  increased  number  of  iterations, 
the  predictor-corrector  solution  was  still  found  to  be  sufficiently  fast 
that  its  use  and  accuracy  (as  seen  from  Tables  1  and  2)  are  practical. 

For  W  =  1.0,  the  successive  approximations  solution  for  the  source  func¬ 
tion  was  found  to  become  very  slow,  with  convergence  not  always  attain¬ 
able  even  after  50  iterations  for  tq  =  5.0,  and  the  accuracy  was  not  good 
when  convergence  was  attained  for  td  =  0.5.  In  general  it  was  found  that 
the  successive  approximations  solution  of  the  integral  equation  became 
very  slow  as  W  -*■  1.0,  and  fine  At  grid  (as  compared  to  the  predictor- 
corrector  method)  was  required  to  achieve  good  accuracy.  Results  for 
W  <  0.  7  by  this  method  were  found  to  converge  in  a  reasonable  amount 
of  computer  time.  The  amount  of  computer  time  required  for  conver¬ 
gence  by  both  the  predictor- corrector  and  successive  approximations 
methods  was  found  to  increase  as  W  increased  and  as  tq  increased.  The 
convergence  speed  should  be  expected  to  decrease  as  tq  increases, 
since  the  Atq  intervals  are  larger.  For  very  large  t0(tq>7)  values  it  is 
necessary  to  use  more  integration  steps  than  the  50  used  in  the 
predictor-corrector  method.  The  Chandrasekhar  method,  on  the  other 
hand,  can  be  shown  to  be  applicable  even  for  a  semi- infinite  medium 
<tq  -*  ®).  For  the  semi- infinite  medium,  p/2  (or  half)  of  the  integration 
constants  in  Eq.  (15)  and  (16)  become  zero. 

Now  the  Chandrasekhar  technique  will  be  utilized  to  compare  the 
results  obtained  by  using  either  the  single  or  double  Gaussian  quadra¬ 
ture  and  also  to  compare  the  agreement  obtained  by  using  different  or¬ 
ders  of  each  type  of  quadrature.  Table  3  shows  the  angular  variation 
of  Phd^^  for  To  =  0-25  and  ri2  =  1.2  and  for  a  range  of  W  for  different 
substrates.  The  results  shown  in  Table  3  were  obtained  by  using 
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Table  3.  Comparison  of  PhdfM1 )  Results  for  the  Single  and 
Double  Gaussian  Quadratures  for  Various  Orders  of 
Quadrature:  n2  =  1.2,  r0  =  0.25 


Quadrature 

Type3 

Quadrature 

Order 

G1 

Aluminum 

W a  0. 4 

Steel 

W  °  0. 7 

Black  Paint 
W- 1.0 

S 

96 

1.714 

0.5442 

0.4127 

0.0857 

S 

80 

2.054 

0.5441 

0.4130 

0.0859 

S 

48 

3.410 

0.5445 

0  4132 

0.0862 

S 

96 

3.934 

0.5442 

04127 

0.0857 

S 

80 

4.716 

0.5443 

0.4128 

00858 

D 

48 

4.773 

0. 5436 

04118 

0.0852 

S 

32 

5.089 

0.5448 

04137 

0.0868 

D 

40 

5.704 

0.5424 

0.4102 

0.0840 

S 

96 

6.169 

0.5438 

0. 4125 

00859 

S 

24 

6.753 

ft  5401 

0.4071 

0.0814 

D 

32 

7.087 

0.5404 

0.4076 

0.0819 

S 

80 

7.397 

0. 5435 

0.4124 

0.0861 

S 

48 

7.833 

0.5436 

0.4128 

0  0865 

s 

20 

8.073 

0. 5419 

0.4102 

0  0844 

s 

96 

8.411 

0.5430 

0.4122 

00862 

s 

16 

10.04 

0.5442 

0. 4145 

0  0887 

s 

80 

10.09 

0.5426 

0.4121 

0  0865 

s 

96 

10.66 

0.5422 

0.4118 

0.0864 

D 

48 

10.96 

0.5415 

0.4109 

0.0858 

D 

20 

11. 14 

0.5448 

0.4172 

0  0910 

S 

32 

11.70 

0.5425 

04127 

0  0875 

S 

48 

12.30 

a  5418 

0.4120 

00871 

S 

80 

12. 79 

0.5412 

0.4115 

00868 

S 

96 

12.91 

0.5411 

0.4114 

00868 

D 

40 

13. 11 

0.5424 

0  4089 

0.0848 

S 

10 

13.12 

0.5393 

04105 

0.0873 

S 

% 

15. 17 

a  5398 

0.4108 

00872 

s 

80 

15.49 

a  5396 

0.4108 

0.0874 

s 

24 

15.55 

0.5359 

04051 

00826 

D 

32 

16.29 

0. 5357 

0.4054 

00832 

S 

48 

16. 79 

a  5390 

0  4108 

0.0880 

D 

48 

17.21 

a  5378 

04093 

0.0870 

S 

96 

17.43 

0.5382 

04101 

00877 

s 

80 

18.21 

a  5377 

04099 

00880 

s 

32 

18.41 

a  5382 

0.4108 

00889 

s 

20 

18.61 

0. 5359 

04075 

00863 

s 

96 

19.71 

a  5364 

0.4093 

0.0883 

D 

40 

20.58 

a  5340 

0.4065 

0.0865 

S 

80 

20.95 

0. 5354 

0.4089 

00888 

D 

10 

21.30 

0. 5263 

0  3954 

0  0780 

S 

48 

21.32 

a  5353 

0.4092 

0.0893 

s 

96 

22.00 

0.5343 

0.4084 

0.0890 

s 

16 

23. 21 

0. 5349 

04105 

00918 

D 

48 

23.49 

0.5322 

0.4068 

00889 

s 

80 

23.70 

a  5327 

0.4077 

0.0897 

24 


A  EDC-TR -73-200 


Table  3.  Continued 


Quadrature 

Type3 

Quadrature 

Order 

fll 

Aluminum 

W-0.4 

Steel 

W-0.7 

Black  Paint 
W- 1.0 

S 

96 

2430 

a  5320 

0. 4074 

0.0898 

S 

24 

2452 

a  5280 

0. 4016 

00850 

S 

32 

25.21 

a  5317 

0.4080 

00911 

D 

32 

25.62 

0. 5271 

04016 

0.0859 

D 

20 

25.68 

a  5344 

0.4124 

0  0950 

S 

48 

25.90 

0.5306 

0.4071 

00909 

S 

80 

26.47 

a  5296 

0.4064 

0.0908 

s 

96 

26.61 

a  5294 

04062 

0  0907 

0 

40 

28.13 

0. 5259 

0.4029 

0  0893 

s 

96 

28.94 

0.5294 

0.4050 

0.0917 

s 

80 

29.26 

0.5262 

0.4049 

00920 

s 

20 

29.45 

a  5244 

0  4025 

0.0902 

D 

48 

29.84 

0.5247 

0.4036 

0.0915 

S 

10 

30.11 

0.5145 

0.3999 

0  0966 

S 

48 

30.55 

a  5248 

04046 

0.0931 

S 

96 

31.29 

0.5234 

0  4037 

0.0930 

S 

80 

32.09 

a  5223 

0.4033 

0.0935 

S 

.32 

32.15 

0.5229 

0.4043 

0  0944 

s 

96 

33.67 

0.5199 

0.4022 

0  0944 

s 

24 

33. 74 

0.5159 

0.3962 

0.0893 

s 

80 

3495 

0.5181 

0.4015 

0.0953 

D 

32 

35.10 

0.5141 

0.3959 

0.0907 

S 

48 

35.28 

a  5178 

0.4017 

0.0959 

D 

40 

35.80 

a  5148 

0.3982 

0.0936 

S 

96 

36.07 

0.5162 

0.4007 

0.0959 

D 

48 

36.27 

0.5152 

0.3996 

0.0954 

S 

16 

36.94 

a  5166 

0.4029 

0.0994 

S 

80 

37.84 

0.5133 

0.3996 

00974 

s 

96 

38.49 

0.5121 

0.3990 

00978 

s 

32 

39. 30 

0.5115 

0.39% 

0.0996 

s 

48 

40.12 

0.5096 

03985 

0.0997 

D 

20 

40.67 

0.5125 

0.4037 

0.1049 

s 

20 

40.75 

0.5065 

0. 3951 

0.0978 

s 

80 

40.79 

0.5081 

03976 

0.0999 

s 

96 

40. 96 

a  5077 

03973 

0.1000 

D 

48 

42.83 

0.5035 

0.3950 

O1011 

S 

24 

43. 38 

a  4990 

0.3893 

0.0970 

s 

96 

43.46 

0.5031 

0.3955 

0.1026 

D 

40 

43.65 

0.5008 

0.3926 

01004 

S 

80 

43. 79 

a  5025 

03954 

01030 

D 

32 

44  87 

0.4965 

0.3888 

00992 

S 

48 

45. 12 

0.5002 

03949 

01050 

s 

96 

46.01 

a  4980 

03937 

01056 

s 

32 

46.78 

0.4973 

0.3944 

0. 1077 

s 

80 

46.87 

0.4964 

0  3933 

0.1068 
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Table  3.  Concluded 


Quadrature 

Type3 

Quadrature 

Order 

61 

Aluminum 

W-0.4 

Steel 

W=a7 

Black  Paint 
W-  1.0 

S 

10 

47.20 

—  0. 4679 

-a  1430 

S 

96 

48.62 

a  4928 

a  1093 

D 

48 

49.62 

a  4900 

0.3903 

a  lioo 

S 

80 

5a  03 

0.4899 

0.3912 

a  ni6 

D 

10 

5a  07 

0.4797 

0.3756 

-*0.0986 

S 

43 

5a  33 

a  4896 

0.3915 

a  1126 

S 

96 

5L30 

0.4872 

0.3903 

a  1138 

D 

40 

51.84 

0.4841 

0.3871 

a  1123 

S 

16 

51.85 

0.4881 

0.3932 

a  1177 

s 

20 

53.00 

a  4819 

a  38  69 

a  1149 

s 

80 

53.30 

0.4832 

0.3894 

a  1178 

s 

24 

53.80 

a  4776 

0.3823 

a  1129 

s 

96 

5407 

0.4815 

0.3890 

0. 1194 

s 

32 

5482 

a  4809 

a  3900 

0.1223 

D 

32 

55.23 

0.4751 

0.3824 

a  1167 

S 

48 

55.86 

0. 4784 

0.3890 

0.1243 

S 

80 

56.72 

a  4832 

0.3883 

a  1261 

D 

48 

56.78 

0.4754 

a  3871 

0. 1253 

D 

20 

56.81 

0.4808 

a  3950 

0. 1321 

S 

96 

56.94 

a  4759 

0.3882 

0.1266 

S 

96 

59.95 

0.4706 

a  3884 

0.1362 

S 

80 

60.34 

0. 4701 

0.3886 

0.1378 

D 

40 

6a  71 

a  4673 

a  3855 

0.1364 

S 

48 

61.91 

a  4682 

0. 3899 

0.1445 

S 

96 

63.15 

0.4663 

a  3903 

0.1496 

S 

32 

63.92 

0. 4664 

a  3925 

0.1548 

S 

80 

6424 

0.4653 

0. 3917 

0. 1552 

D 

48 

64  71 

0.4641 

0.3911 

0. 1567 

S 

24 

66.05 

0.4595 

a  3873 

a  1593 

S 

96 

66.62 

0.4642 

0.  3957 

0. 1694 

D 

32 

67.11 

0.4598 

0.3903 

0. 1670 

S 

20 

67.82 

0.4625 

0.3959 

a  1756 

S 

80 

6a  57 

0.4649 

0. 4010 

0.1840 

S 

48 

68.91 

a  4656 

0.4027 

a  1874 

S 

96 

70.48 

a  4674 

0.4082 

a  2013 

S 

16 

7a  67 

a  4699 

0.4123 

a  2063 

D 

40 

71.31 

0.4670 

0.4089 

0.2073 

S 

80 

73.68 

a  4774 

0.4270 

a  2401 

D 

48 

7457 

0.4811 

0.4329 

a  2526 

S 

96 

75.06 

a  4847 

0. 4383 

0.2615 

S 

32 

76.14 

0.4930 

0. 4505 

0.2820 

s 

48 

7a  49 

a  5154 

0.4808 

0.3335 

D 

20 

79.19 

a  5283 

0. 4981 

a  3568 

S 

80 

80.78 

0.5489 

0.5231 

0.4002 

s 

96 

81.44 

a  5617 

0.5386 

0.4237 

as  - 

Single  Gaussian  Quadrature 

D  -  Double  Gaussian  Quadrature 
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single  Gaussian  quadratures  of  orders  10,  16,  20,  24,  32,  48,  80,  and  96, 
and  double  Gaussian  quadratures  of  orders,  10,  20,  32,  40,  and  48.  Table 
3  shows  no  very  large  or  significant  differences  from  one  order  of  quad¬ 
rature  to  the  other  or  in  the  double  quadrature.  The  points  which  show 
the  biggest  variation  from  the  rest  of  Table  3  are  the  quadrature  direc¬ 
tions  which  transmitted  through  the  coating-vacuum  interface  into  the 
largest  directions  for  p  =  10.  As  noted  by  the  arrows  in  Table  3, 

=  47.20  and  =  50.07  show  the  largest  variation,  and  both  corre¬ 
spond  to  p  =  10  for  the  single  and  double  Gaussian  quadratures,  respec¬ 
tively.  To  clarify  what  is  meant  by  quadrature  order  in  Table  3,  a  32- 
point  quadrature  (double  or  single)  means  16  positive  and  16  negative 
quadrature  directions;  or,  the  quadrature  order  can  be  considered  as 
the  size  of  the  p  x  p  coefficient  matrix  in  Eq.  (10)  (in  this  case,  32  x 
32).  Figure  3  is  a  graphical  representation  of  some  of  the  results  in 
Table  3  and  is  presented  in  order  to  more  clearly  illustrate  the  agree¬ 
ment  between  the  single  and  double  Gaussian  quadratures  for  n2  -  1.2. 
Table  3  and  Fig.  3  show  that  there  is  essentially  no  advantage  with  re¬ 
spect  to  accuracy  in  using  one  type  of  Gaussian  quadrature  over  the 
other.  However,  from  a  computational  point  of  view,  it  is  recommended 
since  both  are  essentially  the  same,  that  the  single  Gaussian  quadrature 


0  10  20  30  40  50  60  70  80  90 

flj.  deg 

Figure  3.  Comparison  of  phd  (p1 )  results  for  various  orders  of  single 
and  double  Gaussian  quadratures. 
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be  used.  There  are  two  reasons  for  this.  First,  for  a  given  order  of 
quadrature,  the  largest  eigenvalue  is  much  smaller  for  the  single  Gaus¬ 
sian  quadrature  than  for  the  double  Gaussian  quadrature,  as  is  shown 
in  Fig.  4.  (The  results  shown  in  Fig.  4  can  be  considered  valid  for  all 
values  of  W  since  the  size  of  the  maximum  eigenvalue  only  changes  by 
about  ±1  on  the  ordinate  scale  in  going  from  W  =  0.0  to  W  =  1.0.  )  This 
means  that  computations  can  be  made  to  a  much  larger  optical  thick¬ 
ness,  Tq,  with  the  single  quadrature  than  with  the  double  quadrature 
since  this  largest  eigenvalue  occurs  as  a  positive  number  in  the  expo¬ 
nent  of  Eqs.  (15)  and  (16);  the  double  quadrature  will  cause  computer 
overflow  at  a  much  smaller  tq  value.  The  second  advantage  is  that  the 
double  quadrature  directions  correspond  to  larger  angles  than  do  the 
single  quadrature  directions.  This  means  that  more  double  quadrature 


Figure  4.  Effect  of  quadrature  order  on  maximum  eigenvalue. 
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directions  are  greater  than  the  critical  angle,  and  thus  p^Ou*)  is  com¬ 
puted  in  fewer  directions  for  the  double  quadrature  since  the  others  are 
trapped  inside  the  coating.  This  situation  becomes  even  more  evident 
as  the  refractive  index  n2  increases.  The  advantage  of  having  more 
quadrature  directions  transmit  through  the  coating- vacuum  interface 
is  that  the  flux  can  be  more  accurately  evaluated  since  the  intensity  is 
known  in  more  directions.  Figure  5  demonstrates  the  difference  be¬ 
tween  the  number  of  transmitted  directions  for  the  single  and  double 
quadratures  an  n2  increases. 


Figure  5.  Effect  of  refractive  index  on  number  of  transmitted 
quadrature  directions. 


Figure  6  shows  that  for  n2  =  1.4  a  16-point  single  Gaussian  quad¬ 
rature  must  be  used  to  achieve  accuracy.  Results  are  shown  for  both 
r0  =  0.25  and  tq  =  5.0.  The  10-point  Gaussian  quadrature  is  seen  to 
yield  results  that  are  too  low.  Whereas  the  10-point  quadrature  showed 


29 


AEDC-TR-73-200 


0.80 

0.70 

0.60 


~  0.50 

3 

'^3 

0.40 


0.30 

0.20 

0 

0  10  20  30  40  50  60  70  80 

Op  deg 

Figure  6.  Order  of  quadrature  needed  to  achieve  accuracy  for  n2  =  1.4. 
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good  agreement  in  Fig.  3  for  n2  -  1.2  (except  for  large  it  is  seen 
that  for  n2  =  1.4  a  higher  order  quadrature  is  needed  for  accuracy. 

This  gives  rise  to  speculation  that  as  n£  increases  the  order  of  Gaus¬ 
sian  quadrature  needed  for  accuracy  also  increases;  Fig.  7  shows  this 
to  be  true.  The  results  shown  in  Fig.  7  were  obtained  by  increasing 
the  quadrature  order  until  two  successive  higher  order  quadratures 
showed  good  agreement.  The  results  given  in  Fig.  7  should  be  con¬ 
sidered  as  the  minimum  quadrature  order  needed  for  accuracy.  Higher 
order  quadratures  than  those  shown  will  likewise  yield  good  accuracy. 
Although  Fig.  7  is  for  the  single  Gaussian  quadrature,  results  were  ob¬ 
tained  for  the  double  Gaussian  quadrature  which  also  showed  the  same 
dependence  of  accuracy  upon  quadrature  order.  Hence,  Fig.  7  may  be 
considered  valid  for  both  the  double  and  single  Gaussian  quadratures. 
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"2 

Figure  7.  Effect  of  coating  refractive  index  on  quadrature  order 
required  for  accuracy. 


5.0  DISCUSSION  OF  ERROR 


The  three  solution  techniques  and  numerical  formulas  used  in  the 
solution  of  the  transport  equation  were  described  in  the  analysis  sec¬ 
tion  of  this  report.  A  discussion  of  the  error  is  now  presented  so  that 
confidence  in  the  results  presented  may  be  achieved.  First,  the 
Chandrasekhar  method  will  be  discussed,  and  then  the  numerical  solu¬ 
tion  and  successive  approximations  solution  will  be  discussed.  From 
the  discussion  of  the  results  in  Table  3  it  is  concluded  that  the  error 
introduced  in  using  one  type  of  quadrature  as  opposed  to  the  other  or 
different  orders  of  quadrature  is  negligible.  Thus  the  error  discussed 
here  will  pertain  to  using  a  fixed  order  of  quadrature.  For  the 
Chandrasekhar  solution,  the  roots  (X2)  of  Eq.  (13)  were  found  by  an 
iterative  procedure,  and  the  X2  values  were  considered  to  have  been 
determined  when  two  consecutive  calculations  were  within  10"^  (double) 
precision)  of  one  another.  Thus  the  eigenvalues  (±X)  were  accurate  to 
10"®.  Once  the  homogeneous  solution,  Eq.  (15)  or  (16),  was  obtained. 
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it  became  possible  for  error  to  be  introduced  in  the  determination  of 
the  p  integration  constants.  In  order  to  estimate  the  error  introduced 
through  the  Cholesky  method,  which  was  used  to  determine  the  integra¬ 
tion  constants,  the  solutions,  Eqs.  (15)  and  (16),  were  put  back  into  the 
boundary  conditions.  The  right  sides  of  Eqs.  (8)  and  (9)  were  subtracted 
from  the  left  sides  for  each  of  the  p/2  quadrature  directions  correspond¬ 
ing  to  each  equation.  This  was  done  for  both  the  double  and  single  Gaus¬ 
sian  quadratures,  for  W  =  0.4,  0.7,  1.0,  for  the  four  substrates  used  in 
Table  1,  and  for  rQ  =  0.25.  It  was  found  that  the  largest  difference  be¬ 
tween  the  two  sides  of  Eqs.  (8)  and  (9)  occurs  at  the  small  optical  thick¬ 
nesses.  Thus,  Fig.  8  is  shown  for  the  small  optical  thickness,  t0  =  0.25. 


0  10  20  30  40  50  60  70  80  90 


Quadrature  Order 

Figure  8.  Effect  of  quadrature  order  on  maximum  error  for 
Chandrasekhar  method,  t„  =  0.25. 
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Figure  8  shows  the  maximum  error  of  all  the  quadrature  directions  as 
a  function  of  quadrature  order;  again,  this  error  represents  the  maxi¬ 
mum  difference  between  the  right  and  left  sides  of  Eqs.  (8)  and  (9)  for 
the  range  of  conditions  described  above.  It  is  seen  here  that  the  double 
quadrature  yields  better  accuracy  than  the  single  quadrature;  however. 
Fig.  8  shows  the  error  for  both  types  and  all  orders  of  quadrature  to 
be  negligible.  The  conclusion  is  that  use  of  the  Cholesky  method  intro¬ 
duced  very  little  error  into  the  final  solution. 

Regarding  the  numerical  solution,  it  should  be  remembered  that 
all  numerical  formulas  are  approximate  since  they  usually  involve  the 
truncation  of  a  Taylor  series.  The  evaluation  of  error  associated  with 
the  numerical  solution-  is  also  approximate.  The  purpose  of  the  error 
term  is  to  offer  a  means  of  obtaining  a  general  idea  of  the  error  com¬ 
mitted  by  numerical  solution  using  a  particular  step  size. 

The  predictor- corrector  method  is  very  convenient  because  it  pro¬ 
vides  a  simple  expression  for  the  error  per  step  estimate.  Call 
Epr(Tn+i)  the  error  in  computing  ipr(Tn+l)  by  using  Eq.  (27);  call 
Eco(Tn+l)  the  error  in  computing  iCo(Tn+l^  by  using  Ech  (25)  to  correct 
this  predicted  value  of  ipr(Tn-t-l)*  Let  Y(7n_j)  be  the  actual  value  of 
i(rn+1);  then 


Ec<,(rn+l^=  “  ’co(rn+P 

(33) 

^pr(rn+P  -  y(rn-f-P  ‘pr(rn+l^ 

(34) 

Subtracting  the  first  equation  from  the  second  gives 

'co^n+P  -  *pr(rn+])  -  ^pr(rn4-P  “  ^co(rn+l^ 

(35) 

The  expressions  for  the  error  terms  of  the  predictor  and  corrector  are 
given  respectively  by 

■  H>>5  iC5>(X,> 

(36) 

and 

E«<'.+1>  -  -^h5i‘5,<X2> 

(37) 

where  Xj  is  contained  in  the  interval  (Tn-3»Tn+l^  anc*  ^2 

is  contained  in 

the  interval  (Tn_  1>  Tn+l)-  Let  tbe  assumption  be  made  that  h  is  suffi¬ 
ciently  small  that  the  variation  between  i(5)  (X^)  and  i(5)  (X2)  is 
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negligible;  then  a  small  error  is  committed  by  using  i(5)  (x)  as  their 
approximate  common  value.  Hence  subtracting  Eq.  (37)  from  Eq.  (36) 
and  replacing  i(5)  (X^)  and  it 5)  (X2)  with  it 5)  (X)  results  in 

VW  -  E..<W  -  If1*5  i!S>  <38> 

Combining  Eqs.  (35)  and  (38)  yields 

|  l>5  -  «..U  (39) 

and  therefore, 

Eco(rn-l^  =  tiprCr^!)  -  ico(rn+  ,)]/29  (  40) 

Equation  (40)  gives  an  approximate  formula  for  the  error  commited  per 
step  of  numerical  integration  by  the  Milne  predictor- corrector  method. 

For  the  purpose  of  plotting  results  which  would  indicate  the  signi¬ 
ficant  general  trends,  it  was  decided  that  the  total  maximum  error  be¬ 
tween  two  successive  iterations  of  i(TQ, /u)  [i  =  1,  ....  p/ 2]  should  be 
less  than  0.001;  thus  the  error  per  step  should  be  less  than  0.001  di¬ 
vided  by  the  total  number  of  steps.  Since  50  steps  were  used,  the  dif¬ 
ference  between  predictor  and  corrector  at  any  station,  t,  for  each 
quadrature  direction  should  be 


'pr(Tn+l^  - 


‘co(rn+l^  -  29  E 


co^rn+ 1^ 


(41) 


and  the  value  of  Eco(Tn+j)  is  approximated  by 


Eco^rn+l^ 


0.001 
No.  steps 


0.001 

50 


0.00002 


(42) 


Thus,  at  each  integration  station  it  was  required  that  ipr  -  ico  <  0.00058. 
Since  the  maximum  optical  thickness  used  was  tq  =  5.0,  this  corresponds 
to  a  "maximum"  value  of  h  =  0.1.  Since  by  Eq.  (37)  the  error  per  step 
is  on  the  order  of  h^,  the  value  of  Eco(rn+j)  “  0.00001,  which  is  less 
than  the  value  required  by  Eq.  (42).  The  values  presented  correspond 
to  the  maximum  error,  and  again  it  should  be  mentioned  that  these  ex¬ 
pressions  give  only  a  reasonable  estimate  of  the  error  committed.  If 
the  difference  between  predictor  and  corrector  is  not  less  than  0.00058, 
the  corrector  may  be  used  again  and  again  until  two  successive  calcu¬ 
lations  are  within  the  given  tolerance.  It  was  not  necessary  to  use  the 
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corrector  more  than  once  except  for  tq  =  5.0,  where  at  a  few  stations 
it  was  necessary  for  multiple  use.  For  the  results  requiring  a  maxi¬ 
mum  error  of  0.001  in  Kt0,£^)  [i  =  1, .  .  .  ,p/2]  between  two  successive 
iterations,  it  was  only  necessary  to  perform  the  calculations  in  single 
precision,  as  will  be  shown  later.  By  an  analysis  similar  to  that, 
just  presented  it  was  possible  to  specify  a  maximum  error  of  10-5  be¬ 
tween  two  successive  iterations  of  i (T0.M^)  [j2  =  l,...,p/ 2] ,  if  the  cal¬ 
culations  were  performed  in  double  precision. 

So  far  only  the  analysis  for  approximating  the  error  associated 
with  the  numerical  formula  has  been  presented.  There  are  also  round¬ 
off  errors  associated  with  any  numerical  calculation  scheme.  An  ef¬ 
fective  procedure  for  determining  the  magnitude  of  the  round-off  error 
is  to  first  carry  out  a  calculation  in  single  precision  and  then  to  use 
double  precision  to  see  if  any  significant  difference  is  present.  How¬ 
ever,  the  formula  errors  and  round-off  errors  are  usually  coupled  to¬ 
gether  in  such  a  way  that  they  cannot  be  separately  approximated. 

Probably  the  most  effective  method  of. approximation  of  the  coupled 
error  is  not  only  to  use  double  and  single  precision  but  also  to  simul¬ 
taneously  reduce  the  step  size  to  one  half  its  previous  value  and  see  if 
any  significant  changes  occur.  In  an  attempt  to  make  sure  that  the  ac¬ 
curacy  of  the  solution  was  within  the  stated  tolerance,  the  numerical 
solution  using  the  10-point  single  Gaussian  quadrature  was  performed 
in  both  double  and  single  precision  and  also  for  various  step  sizes. 

Table  4  illustrates  a  comparison  of  the  results  for  various  step  sizes 
and  for  the  single  and  double  precision  calculations.  As  should  be  ex¬ 
pected,  the  results  at  the  small  optical  thickness  are  much  more  ac¬ 
curate  than  those  of  the  large  optical  thickness  since  the  corresponding 
step  size  is  one-tenth  that  of  the  larger  optical-thickness  step  size.  As 
shown  earlier  in  Table  1,  these  results  are  also  in  good  agreement  with 
the  results  obtained  by  the  Chandrasekhar  solution. 

Also  shown  in  Table  4  is  the  sensitivity  of  the  successive  approxi¬ 
mations  solution  to  the  integration  step  size.  The  results  shown  are 
for  demonstration  purposes,  and  only  single  precision  .calculations 
were  performed.  As  in  the  predictor-corrector  discussion,  the  error 
analysis  was  performed  indirectly  by  varying  the  integration  step  size. 
Hence  the  solution  of  the  integral  equation  is  seen  to  be  essentially  in¬ 
sensitive  to  the  integration  step  size. 
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Table  4.  phd  (ju1 )  Results  for  Various  Step  Sizes  for  the 
Predictor-Corrector  Technique,  n2  =  1.2,  W  =  1 
and  for  the  Source  Function  Method,  n2  =  1.2, 

Milne 

.0, 

W  =  0.4 

Substrate 

Mi 

T0 

Predictor-Corrector.  'A  • 

1.0 

Source  Function.  IV  •  0.4 

Ooubie  Precision 

Sinqle  Precision 

Sinqle  Precision 

50 

Steps 

loo 

Steps 

50 

Steps 

100 

Steps 

200 

Steps 

50 

Steps 

too 

Steps 

Black  Paint 

a  96220 

Q5 

a  1591 

a  1591 

a  1591 

01591 

01591 

004774 

0  04771 

0.79850 

a  1742 

a  1742 

01742 

01742 

01742 

004895 

004892 

a  47403 

a  2294 

a  2294 

02294 

0  2294 

02294 

0  07732 

007729 

n  •  1.48  -  id  00 

0.96220 

5.0 

0.6863 

a  6961 

06878 

0  6878 

0  6877 

007019 

0.06692 

a  79850 

a  7019 

0.7018 

07052 

0.7033 

0.7032 

0  07161 

0.06834 

0.47403 

0.7356 

a  7356 

07367 

0.7368 

0.7366 

0  09893 

0.09576 

Steel 

0.96220 

0.5 

a  5015 

a  5015 

05015 

0  5015 

05015 

02567 

0  2567 

0.79850 

0.5029 

0.5029 

05029 

05029 

osoes 

02348 

0.2347 

0.47403 

0.5192 

05192 

05192 

05192 

0.5192 

02107 

02106 

n  ■  2. 48  -  i3. 43 

Cl  96220 

5.0 

0.7404 

07407 

0  7423 

07423 

07422 

007022 

006695 

0.79850 

0.7533 

a  7535 

0  7550 

0  7550 

0.7549 

0  07162 

006836 

a  47403 

a  7813 

0.7826 

07826 

0  7826 

0  7825 

0.09895 

009577 

Aluminum 

0.96220 

0.5 

0.7472 

0.7472 

0.7472 

07472 

07471 

03721 

03720 

0.79850 

a  7448 

0  7448 

07447 

0  7447 

0  7446 

0.3380 

0  3379 

a  47403 

a  7457 

07457 

0.7457 

0  7457 

07456 

0  2865 

02865 

n-  1.44-i5.32 

0.96220 

5.0 

a  8346 

0  8346 

08346 

08346 

08345 

0  07024 

O0669B 

a  79850 

018428 

0  8428 

08428 

08428 

08427 

0  07164 

0  06838 

a  47403 

0.8606 

0  8606 

08606 

08606 

0  8605 

009896 

0.09579 

Copper 

0.96220 

0.5 

0.9619 

09619 

09613 

09613 

09612 

04475 

0.4474 

0.79850 

a  9611 

09611 

09607 

0  9607 

09606 

04063 

0.4062 

0.47403 

a  9601 

09601 

0  9598 

0  9598 

0  9597 

03394 

0.3393 

n-0.82 -113.00 

a  96220 

5.0 

0.9709 

0  9709 

09706 

09706 

0.9705 

0.07025 

0.06699 

0.79850 

0.9723 

0  9723 

0.9721 

0.9720 

0.9720 

007165 

0.06839 

0.47403 

0.9755 

09754 

0  9752 

0. 9752 

09751 

009897 

0.09580 

6.0  SUMMARY  AND  CONCLUSIONS 


Three  solution  techniques  to  the  radiative  transport  equation  have 
been  demonstrated  on  a  sample  problem,  and  a  comparison  of  the  three 
methods  has  been  presented.  In  addition,  a  comparison  of  different 
orders  of  Gaussian  quadrature  was  presented  for  both  the  double  and 
single  Gaussian  quadrature  approximations.  The  Chandrasekhar  solu¬ 
tion  was  found  to  be  extremely  fast;  only  the  computation  of  the  eigen¬ 
values  required  iteration,  and  this  procedure  converged  rapidly  since 
the  bounds  for  each  eigenvalue  were  known.  The  Milne  predictor- 
corrector  solution  also  converged  very  rapidly  and  yielded  results 
which  agreed  extremely  well  with  the  Chandrasekhar  solution  for  all 
four  substrates,  for  all  values  of  W,  for  all  7q(0-*-5.0),  and  for  all 
viewing  angles.  The  successive  approximations  solution  to  the  integral 
equation  formulation  was  also  found  to  converge  rapidly  for  small  values 
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of  W(W  <L0.70).  However,  for  large  values  of  W(0.70  <  W  <_1.0)  the  in¬ 
tegral  equation  solution  was  found  to  converge  very  slowly,  and  even 
when  a  solution  was  obtained,  its  accuracy  was  not  good.  After  the 
three  solution  techniques  were  compared,  the  Chandrasekhar  method 
was  used  to  compare  various  orders  of  Gaussian  quadratures  for  both 
the  single  and  double  approximations.  Comparison  of  the  double  and 
single  approximations  for  a  fixed  order  of  quadrature  showed  very  little 
difference  in  accuracy.  However,  it  was  found  that  as  the  coating  re¬ 
fractive  index  increased,  the  order  of  the  Gaussian  quadrature  had  to 
be  increased  to  maintain  accuracy.  From  a  computational  point  of 
view,  the  single  Gaussian  quadrature  was  seen  to  have  two  advantages 
over  the  double  quadrature.  First,  the  largest  eigenvalue  is  consider¬ 
ably  smaller  for  the  single  quadrature  approximation,  thus  allowing 
computation  to  much  larger  tq  values  before  encountering  computer 
overflow.  Second,  the  double  quadrature  directions  correspond  to 
larger  angles  measured  from  the  substrate  normal.  This  means  that 
more  quadrature  directions  are  trapped  inside  the  coating  past  the  crit¬ 
ical  angle;  the  result  is  that  the  single  quadrature  approximation  has 
more  quadrature  directions  which  transmit  through  the  coating-vacuum 
interface.  Hence  more  information  is  transmitted  in  more  directions. 
This  implies  integration  with  respect  to  /Li1  in  determining  that  flux  can 
more  accurately  be  performed. 

Finally,,  a  discussion  of  error  for  each  technique  was  given,  and 
the  sensitivity  of  the  predictor- corrector  and  successive  approxima¬ 
tions  solutions  to  finite-difference  step  size  was  shown.  All  three  meth¬ 
ods  presented  are  easy  to  understand  and  use;  also,  the  Chandrasekhar 
method  can  be  used  as  an  initial  guess  criterion  for  the  other  techniques 
when  W  is  a  function  of  t.  Furthermore,  all  three  techniques  can  be 
extended  to  include  emission  (particular  solutions)  and  anisotropic 
scattering. 

The  results  discussed  here  are  important  in  choosing  which  tech¬ 
nique  to  use  when  solving  the  radiative  transport  equation.  The  solu¬ 
tion  to  the  transport  equation  is  a  powerful  tool  to  be  employed  in  de¬ 
termining  the  reflectance  and  transmittance  behavior  of  absorbing  and 
scattering  condensed  gas  coatings  as  a  function  of  viewing  angle,  coat¬ 
ing  thickness,  absorption,  scattering,  and  substrate  effects.  Also,  the 
solution  to  the  transport  equation  can  be  used  in  conjunction  with  ex¬ 
perimental  data  for  determining  the  coating  thickness  and  optical  prop¬ 
erties  of  condensed  gases. 
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